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(57) ABSTRACT 

Cartesian mesh generation is accomplished for component 
based geometries, by intersecting components subject to 
mesh generation to extract wetted surfaces with a geometry 
engine using adaptive precision arithmetic in a system which 
automatically breaks ties with respect to geometric degen- 
eracies. During volume mesh generation, intersected surface 
triangulations are received to enable mesh generation with 
cell division of an initially coarse grid. The hexagonal cells 
are resolved, preserving the ability to directionally divide 
cells which are locally well aligned. 

35 Claims, 5 Drawing Sheets 
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TRIANGLE GEOMETRY PROCESSING FOR 
SURFACE MODELING AND CARTESIAN 
GRID GENERATION 

This application claims the benefit of provisional appli- ' 
cation No. 60/068,846, filed Dec. 12, 1997. 

ORIGIN OF INVENTION 

The invention described herein was made by an employee } 
of New York University and employees of the United States 
Government and may be manufactured and used by or for 
the Government for governmental purposes without pay- 
ment of any royalties thereon or therefor. The invention 
described herein was made in the performance of work ( 
under Air Force Office of Scientific Research (AFOSR), 
Department of Energy (DOE), and NASA contracts and is 
subject to the provision of Section 305 of the National 
Aeronautics and Space Act of 1958. Public Law 85-568 (72 
Stat. 435; 42 U.S.C. Section 2457). : 

FIELD OF THE INVENTION 

This invention relates to computer technology and. in 
particular, to computer applications for processing of trian- 
gulated surfaces and Cartesian mesh generation. 

BACKGROUND OF THE INVENTION 

Adaptive Cartesian mesh solutions have been attempted 
for problems involving complex geometries. Flow solvers 
and mesh generation schemes for use with arbitrary' geom- 
etries have been sought. Some approaches have included 
beginning with a surface triangulation constrained to the 
intersection curves of particular components. In contrast, a 
component based approach requires that each element of a 
geometry be described as a single closed entity. However, 
using the component based approach adds complexity dur- 
ing grid generation. Since components may overlap, Carte- 
sian cell surface intersections must be characterized and 
those internal to components must be rejected. Such rejec- 
tion filtering increases complexity. Assembling surface 
geometries and creating volume grids on a computer (for 
example, for computational physics simulations) is accord- 
ingly a difficult and time-consuming process. Thus, there is 
a need to improve accuracy and efficiency of the comput- 
erized surface creation, component definition, and volume 
mesh generation. It is desirable that an automatic method for 
discretizing the domain around arbitrarily complex geom- 
etries be accomplished. 

SUMMARY OF THE INVENTION 

According to the present invention, computer techniques 
for processing of triangulated surface geometry and Carte- 
sian mesh generation is made particularly robust and effi- 


2 

those from computational geometry and computer graphics, 
are used to generate a body-intersecting Cartesian volume 
mesh that discretizes the surrounding flow field. 

According to the present invention, the creation of com- 
plex surface definitions is simplified. Further according to 
the present invention, selected complex surface definitions 
are made accurate and precise to a degree not known in the 
prior art. Accordingly, accurate computational aerodynamic 
performance estimation is facilitated. Additionally, Carte- 
sian grid generation is made more effective. More 
° specifically, the computer system according to the present 
invention forms the Boolean union of multiple-component 
geometries, and eliminates selected portions of the 
assembled geometry that are w ithin the assembled geometry 
and thus hidden from exterior view. According to the present 
s invention, computations of surface-related quantities are 
computed w'ith increased accuracy, enabling the effective 
production of data images describing features such as lor 
example, w ithout limitation, the exposed sui t ace area and 
the internal volume of a selected structure, to produce a 
n structure descriptive information dataset representing pro- 
cessed surface data. The processed surface data is then 
prov ided to the volume mesh generator for further compu- 
tation. A surrounding flowlield for the selected structure is 
automatically discretized into rectangular hexahedra (cells), 
with the size of the hexahedra being reduced in areas of large 
surface curvature or strong flow’ field gradients. Cells that 
are split into multiple distinct and unconnected regions by 
thin pieces of the surface are identified, and the individual 
regions are marked and saved in the computer memory for 
use, for example, in flow field simulations. According to the 
present invention, a Cartesian mesh generation system is 
provided that is capable of identifying and correcting cells 
that are split into distinct, unconnected regions. For example 
and without limitation, the disclosed technology may be 
employed for applications that include aerodynamic perfor- 
mance estimation and optimization for ground and aero- 
space vehicles; computer-based geometry modeling (CAD/ 
CAM) systems; rapid prototyping systems 
(sterolighography); computer-based visualization systems 
(medical, chemical, and other imaging); computational 
physics modeling systems (CEM/CFD, etc.); semiconductor 
device modeling; and internet and data transfer applications 
(including substantial reductions in the size of VRML and 
other geometry-related datasets). 

BRIEF DESCRIPTION OF THE DRAWINGS 
FIG. 1 is a diagram of a refinement tree according to one 
embodiment of the present invention; 

FIG. 2 a is a diagram of a method according to the present 
invention for dividing a single cell of a selected structure 
into multiple polyhedra by intersected geometry techniques; 
FIG. 2b is a diagram of a method for polygon matching; 
FIG. 3A is a flow chart of a method of intersection 
determination according to one embodiment of the present 


cient. In particular according to the present invention, com- 
plex two and three dimensional computer configurations are 
produced from combinations of simpler components. Por- 
tions of components that lie inside of other components are 
automatically trimmed away according to the present 
invention, leaving only the external (exposed) portions of 
the selected component collection. The intersections 
between components are accurately expressed in the result- 


invention; 

FIG. 3B is a flow chart of a method of intersection 
determination according to another embodiment of the 
present invention; 

FIG. 4 is a flow chart of a method of grid generation 
according to one embodiment of the present invention; and 

FIG. 5 is a diagram of nomenclature of a non-uniform 
mesh in 1-D. 


ant triangulation. 

According to one embodiment of the present invention, 
geometric degeneracies (tie breaking) are unambiguously 65 
resolved. Efficient data structures are selectively stored in 
computer memory, and specialized techniques, including 


DETAILED DESCRIPTION OF A PREFERRED 
MODE 

FIG. 1 is a diagram of a refinement tree according to one 
embodiment of the present invention. The disclosed embodi- 
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merits are implemented by software running preferably on new- children are tagged for refinement in both Y and Z, and 
general purpose computer hardware. For example, the dis- they accordingly have codes (0 1 1 ). After each of these cells 
closed system may employ hardware available in typical are divided next in Y, there will be 4 children, and each of 
workstations and servers, such as an SGI (Indigo, Onyx and these inherits a code of (0 0 I ), implying that there is still a 
0iigin_000) or SUN ( Sparc/Ultraspare ) desktop machine. 5 division in Z remaining. After division in Z is complete. 
The present invention may also be implemented on a com- there will be 8 children cells, and each will have a code of 
putei such as Cray C90. Alternatively, the system and (() 0 0) indicating that refinement is complete. The refine- 
method ot the present invention may also be implemented to ment tree according to the present invention is illustrated in 
run on a Pentium PC hardware system. According to one FIG. 1. Each refinement according to the present invention 
embodiment of the present invention, bitmasks are used for 10 constitutes a simple binary partifion of the cell. Thus, the 
uniquely locating/constraining points/faces/edge/comers on refinement is implemented with a single cell division algo- 
a cell. In particular, three dimensional points associated with rithm which splits a cell into two children. According to the 
each Cartesian cell are assigned a boolean code with a length present invention, split cells of adjacent polyhedra are 
equal to _ times the number of dimensions ot the cell (in 3D. matched by testing Cartesian edge fragments for overlap. A 
there are 6 bits, similar to the OUTCODE). This bitmask is 15 single Cartesian cell is accordingly divided into several 
used to uniquely constrain the location of the 3D points to polyhedra by intersected geometry (see cell A in FIG. 2a). 
certain parts of the Cartesian cell. Points within the cell are When this occurs, a face (which is a polygon) on the original 
assigned a Zero. The six positions for this bitmask for 3D “facepoly AT* corresponds to a “facepoly A2" in FIG* 2a. 
Caitesian cells are: The face polygons one each cell accordingly have corre- 

(Lovv_X. High_X, Lovv_Y, High_Y. Low_Z, High__Z) 20 sponding polygons on the neighboring cell (cell B in FIG. 

It a 3D point is “ON’* (constrained to) a face of the cell, 2a). However, when a cell A is split into multiple polyhedra 
the corresponding bit is set to I otherwise it is set to “O' . by the intersected geometry, the correspondence of its face 

Points ON edges of the cell may be thought of as being polygons with those of its neighboring cells mav be lost. 

“ON two laces of the cell, and therefore, they will have 2 This is a specific case of “polygon matching." According to 

bits set. For example, if a point has its Low_X and Low__Y 25 this process, two-dimensional polygons (facepolys like Al, 
bit set, then that point must lie precisely ON the edge formed A2. B 1. B2 are planar and therefore may be constructed in 
by the Low__X and Low__Y faces of the Cartesian cell. two spatial dimensions only) are matched with correspond- 
Similarly points constrained to the comers of the Cartesian ing two dimensional polygons on abutting polyhedra in three 

cell have 3 bits set. since at cell comers, three faces share a dimensional space. Polygon matching is performed accord- 

common point. 3D points interior to the cell have no bits set. 30 ing to the present invention during Cartesian mesh genera- 
Therefore, the number of bits set in this mask uniquely tion in operations not only for processing split cells, but also 
distinguish points that lie on edges, on faces, or on comers during mesh refinement boundary processing. Thus, accord- 

of the Caitesian cell. The positions of the bits which are set ing to the present invention, an intersection test for two 

further classify the location of the 3D points with respect to planar polygons which exist in a three dimensional domain 

the cell. This system of classification makes it possible to 35 is performed. Floating point inaccuracy makes typically 
rapidly compare locations of objects, and to do constrained employed algorithms problematic to implement robustly for 

comparisons, such as intersection testing of two line seg- such operations. Furthermore, all of the points which lie on 

ments which are constrained to a Cartesian cell face. This Cartesian edges or lie on the faces of the Cartesian cell are 

classification frequently permits intersection testing to be “constructed" geometry, and are therefore subject to round- 

earned out in lower dimensions. According to one embodi- 40 off error from their construction using floating-point math, 
ment of the present invention, bitwise tags are used for Such difficulties eliminate the approach of trying to match 

directional refinement. In particular, Cartesian cells are points within one polygon with points within the other, or 

refined or divided in any combination of Cartesian direc- casting a ray from within one polygon and seeing if it 

tions. For example, a Cartesian cell may be divided in one intersects the boundary of the other. The approach according 

direction (i.e., X, Y, or Z only), which would partition the 45 to the present invention thus simplifies into a matching test 
original parent cell into 2 children. In other instances, a cell in a single dimension. Polygons are therefore processed by 
may require refinement into two directions ( i.e. X&Y, X&Z, matching segments of their edges which are coincident with 
or Y&Z). Division into two directions splits a single parent edges of the Cartesian cell of interest. The bitmasks 
cell into 4 children. Finally, a cell may be divided in all three employed for the vertices of the polygon being processed 
directions, X, Y, and Z. Such a cell would produce 8 children so make it possible to unambiguously identify which polygon 
cells after division. To specify the division required accord- edges are coincident with particular edges of the Cartesian 
ing to cell, a 3 bit code is used, with each bit corresponding cell. When the bitmasks indicate that an edge of a polygon 
to division in a Cartesian direction. Any bit ordering con- A is constrained to the same Cartesian edge as an edge of a 
vention would work, but take for example (bit_X,bit_Y, polygon B, there is accordingly a possibility that polygon A 
bit— Z). A set bit indicating division in that direction is 55 matches polygon B. Constraining to the Cartesian edge 
accordingly required. The code (10 0) thus signifies that the removes two coordinate directions from the overlap test, and 
cell is to be divided in X only (with two children cells). there is only one dimension left to check for a match. 
Similarly, the code (10 1 ) would indicate division is Therefore, matching the face polygons of one cell with the 
required in X and Z directions ( with 4 children cells). A code corresponding face polygons of another cell reduces to an 
of (1 1 1 1 uniquely specifies a single parent cell. When cell 60 overlap test of two line segments which are constrained to an 
division occurs according to the present invention, it occurs edge of a Cartesian cell. Furthermore, this 1-D test is robust, 
on a direction-by-direction basis. At the time of division, the since non-matching polygons will occupy no overlapping 
bit for the current division gets unset (assigned a value of segments of the Cartesian edge. 

“0"). After a cell is divided, its code is inherited by its FIG. 2B illustrates an algorithm according to one embodi- 
children. Thus, if a cell is initially tagged with a code (11 65 ment of the present invention, for accomplishing polygon 
1 ), it is divided in X. This unsets the X bit, and the remainder matching. In particular, the bitmasks of a polygon_a and a 
of the code is inherited by the two children. Each of these polygon_b determine that edge AB of polygon a and the 
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edge CD of polygon_b lie along the same ray. Furthermore, 
in the 2D plane defined by the polygons, both polygon__a 
and polygon_b lie above that ray. Thus, if the edge AD 
overlaps with the edge CD, the polygons intersect. This 
procedure is also used for matching faces of surface inter- 
secting cells at refinement boundaries, and for split-cells at 
refinement boundaries. 

FIG. 3A is a flow chart of a method of component 
intersection determination 29 according to one embodiment 1( ) 
of the present invention. In particular, component intersec- 
tion determination 29 includes the processing of triangulated 
surface geometries by locating 31 triangle intersection can- 
didates with an alternating digital tree (ADT) data structure 
or use of another spatial data structure such as for example 1:1 
a binary tree, a K-d tree, a regional quad/octree, a bucket 
quad/octree, a range tree, a Kdb-tree, an R-tree, an hB tree, 
an MX tree, a bin tree, a sector tree, a 2n-tree. sequential and 
inverted lists, and linear trees as well as other techniques as ->o 
described in Samet, H., The Design ami Analysis of Spatial 
Data Structures, Addison -Wes ley Series on Computer Sci- 
ence and Information Processing: Addison- Wesley Publish- 
ing Company. (1990), according to the present invention. 

Use of the ADT data structure enables implementation of a 25 
proximity query according to the present invention, in which 
a list of candidate triangles is reduced, saving considerable 
processing time. The ADT is a hyperspace search technique 
which converts the problem of searching for finite sized 
objects in d dimensions to the simpler one of partitioning a 30 
space with 2d dimensions. Since the searches are not con- 
ducted in physical space, objects which are physically close 
together are not necessarily close in the search space. 
Accordingly, to improve lookup times, a component bound- 
ing box filter is applied to the triangles before inserting them 
into the tree. Triangles which are not in the bounding box of 
a component other than its ow n, are not inserted into the tree. 
This filtering technique not only reduces the tree size, but 
also improves the probability of encountering an intersection 4( 
candidate within the tree, since the ADT data structure is not 
crowded with irrelevant geometry'. The component intersec- 
tions which are actually checked according to the present 
invention are only those relating to candidate triangles from 
the list of spanning triangles returned from the tree. With the 4f 
task of intersecting a particular triangle reduced to an 
intersection test between that triangle and those on the list of 
candidates provided by the ADT, the intersection determi- 
nation according to the present invention reduces to a series 
of triangle to triangle intersection determinations. According 5( 
to the present invention, each intersecting triangle to triangle 
pair contributes a single segment to the final polyhedra that 
will comprise the wetted surface of the particular configu- 
ration being evaluated. Triangle to triangle intersections 
according to the present invention are determined by a ‘ 
Boolean logic check. For two triangles to intersect in three 
dimensional space, two edges of one triangle must cross the 
plane of the other triangle, and there must be two pierced 
triangle edges. The Boolean logic check according to the ^ 
present invention includes determination of a triple product 
without division. More particularly, the Boolean primitive 
for the 3D intersection of an edge and a triangle is based 
upon the signed volume of a tetrahedron in R 3 which is 
defined by a segment intersecting a triangle, which can be 6 
determined by evaluating the determinant of the following 
matrix: 


</<> a | (i; I 

h n h x hi 1 
n, r, <■: 1 
d u <l\ </: I 

where a, b, c, and d are the vertices of the determinant, and 
aO, a 1 , and a2 are the coordinates of the vertex a; bO. b 1 . and 
b2 are the coordinates of the vertex b; cO, cO. and c2 are the 
coordinates of the vertex c; and dO, dl. and d2 are the 
coordinates of the vertex d of the tetrahedron being evalu- 
ated. The determinant is positive w hen points a, b, and c are 
arranged in a counterclockwise circuit on a plane viewed 
from a point d. Positive and negative volumes define the two 
states of the Boolean test, while zero indicates that the points 
a. b, c, and d are coplanar. If the vertices are coplanar. this 
constitutes a “tie" which is resolved with a tie-breaking 
algorithm according to the present invention. A linked list is 
produced of all segments identified which result from inter- 
sections between triangle to triangle pairs. The pierce points 
at which the segments intersect the triangles are then deter- 
mined conventionally. The intersecting segments divide the 
intersected triangles into polygonal regions which are either 
inside or outside the component body subject to operation. 
The portions of the polygon which are inside the body are 
triangulated, and the triangles which lie inside the body are 
rejected. According to one embodiment of the present 
invention, a well known two dimensional Delaunay algo- 
rithm is used within each original intersected triangle. Next, 
the results of the exact arithmetic are filtered 32 using 
geometric primitives to determine whether a particular tri- 
angle is internal to the wetted surface of a body configura- 
tion being evaluated and is to be deleted. Deletion of 
triangles which are internal to the configuration is accom- 
plished according to the present invention, by ray casting. It 
is determined whether a point is within a particular compo- 
nent of a selected configuration if a ray cast from the point 
intersects the closed boundary of the configuration an odd 
number of times. Ray casting according to the present 
invention is accomplished with a Boolean operation. Count- 
ing the number of intersections determines a triangle's status 
as inside or outside. According to the present invention, the 
status of a tested triangle as inside or outside is passed on, 
i.e., painted, onto each of the three adjacent triangles sharing 
its edges. This passing of status propagates until a con- 
strained edge is reached. According to one embodiment of 
the present invention, the painting is implemented with a 
stack, to avoid the overhead associated with recursion. Then, 
degeneracies are identified and resolved 33 . A degeneracy is 
defined as an outcome from determinant evaluation which 
exactly equals zero, in connection with ray casting degen- 
eracies and during Boolean triangle-to-triangle intersection 
determination. Once a degeneracy has been determined, an 
adaptive precision, floating point filtering and virtual per- 
turbation technique according to one embodiment of the 
present invention is implemented. Accordingly, the deter- 
minant evaluation cases are established to segregate cases in 
which exact determinant evaluation is unnecessary, because 
the determined evaluation calculation produces a result 
which is greater than the error bound established for deter- 
minant evaluation. If the determinant evaluation results in a 
value which is less than the error bound, then the result 
cannot be relied upon, and the exact arithmetic evaluation 
calculation needs to be undertaken. The identification of 
degeneracies, in which its is unclear whether a point, i.e., a 
triangle vertex, is within or outside of a configuration, is 
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determined according to the present invention, by evaluating 
the determinant representing the volume of a tetrahedron 
defined by a point off a plane surface holding a triangle 
intersected by a segment from the point being evaluated, 
with a floating point calculation and comparing the results of 5 
the calculation with an error bound for the floating point 
calculation, as per an adaptive precision exact arithmetic 
procedure developed by Shewchuk. 1 R. See his article 
“Adaptive Precision Floating-Point Arithmetic and Fast 
Robust Geometric Predicates.'* CMU-CS-%-140. School of to 
Computer Science, Carnegie Mellon Univ., 19%. If the error 
bound exceeds the volume determined by floating point 
calculation, then the result is degenerate, meaning that the 
point being evaluated is substantially present on the surface 
of the configuration subject to evaluation. When a degener- is 
ate geometry is identified by the exact arithmetic approach 
just indicated, a simulation of simplicity approach is under- 
taken to apply a virtual perturbation to break a tie which 
placed the point under evaluation substantially on the sur- 
face of the configuration being evaluated. Accordingly, final 20 
triangulated wetted surface regions are identified 34. and ray 
casting is minimized by restricting ray casting to a prede- 
termined bound or region. 

FIG. 3B is a flow chart of a method of intersection 
determination according to another embodiment of the 25 
present invention. In particular. FIG. 3B is a flow chart of a 
method of component intersection determination 39 accord- 
ing to one embodiment of the present invention. The com- 
ponent extraction method according to the present invention 
permits robust extraction of a resultant wetted surface from 30 
a plurality of intersected components. In particular, compo- 
nent intersection determination 39 includes the filter 
processing, i.e., filtering, of triangulated surface geometries 
by performing triangle intersection 41 subject to use of an 
alternating digital tree ( ADT) data structure, according to the 35 
present invention. Performance of triangle intersection 41 
includes determination of which triangles are candidates for 
intersection. Each component surface has been pre- 
triangulated, and the intersection of multiple pre- 
triangulated component surfaces produces a single wetted 40 
surface comprising subsets of triangles from the respective 
original component surfaces and additionally new triangles 
from triangulation operations, according to the present 
invention. ADT filtering of triangles permits identification of 
which triangles of the original component surfaces are 45 
relevant to the establishment of new triangulations (i.e., 
retriangulations of intersected triangles of intersecting com- 
ponent surfaces) for the combined wetted surface. The 
triangles are filtered according to one embodiment of the 
present invention with a ADT data structure which enables so 
implementation of a proximity query according to the 
present invention. Accordingly, the triangle filtering permits 
the combined list of candidate triangles from intersected 
components to be reduced, saving considerable processing 
time in the determination of which component triangles are 55 
actually intersecting. The ADT method is a hyperspace 
search technique which converts the problem of searching 
for finite sized objects in d dimensions to the simpler one of 
partitioning a space with 2d dimensions. Since the searches 
are not conducted in physical space, objects which are 60 
physically close together are not necessarily close in the 
search space. Accordingly, to improve lookup times, a 
component bounding box filter according to the present 
invention is applied to the triangles before inserting them 
into the tree for analysis. Triangles which are not within an 65 
established bounding box of a component other than its ow n, 
are not inserted into the tree for further processing. This 
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filtering technique not only reduces the tree size, but also 
improves the probability of encountering an intersection 
candidate within the tree, since the ADT data structure is not 
crowded with irrelevant geometry' elements. The component 
intersections which are actually checked for intersection of 
triangles, according to the present invention, are only those 
candidate triangles from the list of spanning triangles, which 
are returned from the tree analysis. Thus, the task of inter- 
secting a particular triangle is reduced to an intersection test 
between that triangle and those other triangles on the list of 
candidates provided by the ADT, and the intersection deter- 
mination according to the present invention reduces to a 
series of triangle-to-triangle intersection determinations. 
According to the present invention, each intersecting 
triangle-to-triangle pair contributes a single segment to the 
final polyhedra that will comprise the ultimate desired 
wetted surface of the particular configuration being evalu- 
ated. Triangle-to-triangle intersections according to the 
present invention are determined by a Boolean logic check. 
For two triangles to intersect in three dimensional space, two 
edges of one triangle must cross the plane of the other 
triangle, and there must be two pierced triangle edges. The 
Boolean logic check according to the present invention 
includes determination of a triple product without division. 
More particularly, the Boolean primitive for the 3D inter- 
section of an edge and a triangle is based upon the signed 
volume of a tetrahedron in R * which is defined by a segment 
intersecting a triangle, which can be determined by evalu- 
ating the determinant of the following matrix: 

Oil a I U-y I 

bn b\ b 2 1 

to c i c 2 1 

do cl | d 2 I 

where a, b, c, and d are the vertices of the determinant, and 
aO, al, and a2 are the coordinates of the vertex a; bO, bl, and 
b2 are the coordinates of the vertex b; cO, c 1 , and c2 are the 
coordinates of the vertex c; and dO, dl, and d2 are the 
coordinates of the vertex d of the tetrahedron being evalu- 
ated. Next, the results of the exact arithmetic are filtered 32 
using geometric primitives to determine whether a particular 
triangle is internal to the wetted surface of a body configu- 
ration being evaluated and is to be deleted. Deletion of 
triangles which are internal to the configuration is accom- 
plished according to the present invention, by ray casting. It 
is determined whether a point is within a particular compo- 
nent of a selected configuration if a ray cast from the point 
intersects the closed boundary of the configuration an odd 
number of times. Ray casting according to the present 
invention is accomplished with a Boolean operation. Count- 
ing the number of intersections determines a triangle's status 
as inside or outside. According to the present invention, the 
status of a tested triangle as inside or outside is passed on, 
i.e., painted, onto each of the three adjacent triangles sharing 
its edges. This passing of status propagates until a con- 
strained edge is reached. According to one embodiment of 
the present invention, the painting is implemented with a 
stack, to avoid the overhead associated with recursion. Then, 
degeneracies are identified and resolved 33. The identifica- 
tion of degeneracies, in which it is unclear whether a point, 
i.e., a triangle vertex, is within or outside of a configuration, 
is determined according to the present invention, by evalu- 
ating the determinant representing the volume of a tetrahe- 
dron defined by a point off a plane surface holding a triangle 
intersected by a segment from the point being evaluated. 
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with a floating point calculation and comparing the results of combined with AND and OR binary primitives, to determine 

the calculation with an error bound for the floating point whether a segment spanning the vertices represented by the 

calculation, as per an adaptive precision exact arithmetic respective outcodes intersects a cell boundary. Questionab e 

procedure developed by Shewchuk. J. R. See his article bits are screened out before evaluation, and are not used or 

“Adaptive Precision Floating-Point Arithmetic and Fast 5 determining intersections. Then, split cells in adjacent poly- 
Robust Geometric Predicates." CMU-CS-%-140, School of hedra are matched 145 by testing adjacent Cartesian edge 

Computer Science. Carnegie Mellon Univ.. 1 996. If the error fragments for overlap. Thereafter, the polyhedra in split cells 

bound exceeds the volume determined by floating point are identified 146 , employing tags 147 tor directional refine- 

calculation, then the result is degenerate, meaning that the ment. 

point being evaluated is substantially present on the surface to FIG. 5 is a diagram of a non-unitorm, one-dimensiona 
of the configuration subject to evaluation. When a degener- mesh. In particular, there is shown an irregu ai one- 
ate ^eometrv is identified by the exact arithmetic approach dimensional mesh which can be represented mathematically 
iustlndicated, a simulation of simplicity approach is under- according to well-known techniques. The present invention 
taken to apply a virtual perturbation to break a tie which applies to single and multiple dimensional meshes and mesh 
placed the point under evaluation substantially on the sur- 15 representations. With such mesh schemes, it is possible to 
face of the configuration being evaluated. Accordingly, final produce diagrams of a Cartesian mesh for an attack heli- 

trianaulated wetted surface regions are identified 34 , and ray copter configuration and a close-up ot the mesh through the 

casting is minimized bv restricting rav casting to a prede- left wing and stores. The present invention facilitates pro- 

termined bound or region. " Auction of Cartesian mesh representations of structural con- 

FIG 4 is a flow chart of a method of Cartesian volume 20 figurations made of multiple components such as or 
mesh generation 139 . beginning with final wetted surface example, without limitations, an attack helicopter and its 
trian filiation according to one embodiment of the present structural components including wings and stores. With sue 
invention. In particular? mesh generation 139 includes prox- mesh schemes, it is possible to produce diagrams of isobars 
imitv testing 140 of surface triangles from the final wetted resulting from inviscid flow analysis of an attack ehcopter 
surface triangulation, to enable determination of surface 25 configuration determined on a mesh with 1.2M cells. The 
intersections with particular Cartesian cells, initially, a pre- present invention facilitates production of Cartesian mesh 
determined number of surface mangles are inserted into an representations of such configurations made ot multiple 
alternating digital tree (ATD) or using another spatial data components such as for example, without limitations, an 

structure such as for example a binary tree, a K-d tree, a attack helicopter and its structural components including 

regional quad/octree, a bucket quad/octree, a range tree, a 30 wings and stores. With such mesh schemes, it is possible to 

Kdb-tree, an R-tree, an hB tree, an MX tree, a bin tree, a produce a scatter plot of mesh size vs. computation time tor 

sector tree, a 2”-tree, sequential and inverted lists, and linear a double-teardrop configuration. As shown, computation 

trees as well as other techniques as described in Samet, H., time increases linearly with the number oi cells represented. 

The Design and Analysis of Spatial Data Structures , FIG. 1-5 of Appendix B are diagrams ol an AMR mesh 

Addison-Weslev Series on Computer Science and Informa- 35 blocking for cells tagged around a circular region in - ' 
tion Processing ; Addison-Wesley Publishing Company, respectively with cells tagged for refinement with patch 
( 1990) The triangles are tested for intersection with Carte- outlines indicated by heavy lines, and according to a re er- 

sian cells, and if an intersection is found, the cell is enced sample calculation. In essence, AMR is a hierarchical 

subdivided, with the child cells inheriting the triangle list of inter-mesh communication scheme. It relies upon block - 

its parent. Initially, a cell will test for intersection with the 40 structured, locally refined mesh patches to increase t e 
entire list of triangles. A divided cell will only test for resolution of an underlying structured coarse gnd. Mesh 

intersection with the triangles with which its parent cell had patches at different levels are stored in a tree-based hierar- 

an intersection relationship, according to one embodiment of chy. The method begins with an estimate of the local 

the present invention. The triangle lists connected to a truncation error within each cell of a coarse gnd patch. Cells 

surface intersection TcuO Cartesian cell get shorter by 45 with high error are then tagged for later refinement. Rather 
approximately a factor of 4 with each subdivision. The than divide each individual cell, the tagged cells are orga- 

Cartesian cells are subject to directional division for a nized into rectangular grid patches which usually contain 

predetermined number of times, beginning with an initially hundreds to thousands of cells per patch. Accordingly, it is 

coarse grid of hexagonal cells. To establish whether a possible to produce a diagram of density contours for a time 

Cartesian cell is a candidate for division, the cell is tested for 50 dependent shock reflection from a double ramp using AMR 
intersection with a plurality of points constituting vertices of on Cartesian meshes, with mesh patch outlines being shown, 

triangles from the final wetted surface triangulation provided according to the prior art. In particular, FIG. 6 shows an 

prior to mesh generation. Each vertex is assigned an outcode example of a 2-D time dependent double ramp AMR reso- 

associated with its location with respect to a particular lution. . . 

Cartesian cell, to enable a fast intersection test. Determina- 55 FIG. 1-7 of Appendix B are respective diagrams of a 
tion of rapid triangle Cartesian cell intersections with an parent-child relationship in a Cartesian mesh using a 

integer coordinate system is accomplished 141 by cell regional quadtree data structure, and a corresponding sub- 
division of intersected cells in a multilevel Cartesian grid tree and a scheme for array based storage. The Figure 

142 . Surface intersecting Cartesian cells are automatically particularly shows an example of an ordering for a single 

refined, i.e., divided, a predetermined number of times, and 60 parent node and its four children in a 2-D quadtree mesh, 
the refinement is propagated to neighboring cells by a paint The quadtree shown is referred to as a regional quadtree, 

operation. Further refinement includes curvature detection since each node in the tree corresponds to a region in the 

including detection of the angular variation of a surface domain. Each node in the tree may be pointed to by only one 

normal. In particular, bitmasks are used 143 for uniquely parent node, and only the root has no parent. Leai nodes are 

locating cell features, and outcode data structures are used 65 distinguished by the fact that they have no children. The 
144 for rapid intersection determination between selected Figure additionally shows a simple 1-D array-based storage 

triangles and cell segments. The outcodes are logically scheme for mapping a regional quadtree tree to memory. 
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FIG. 1-8 of Appendix B are respective diagrams of a tree of the same number of objects. As a resuit, fewer 
quadtree mesh, a regional quadtree structure and an array operations are required when the tree is used to retrieve data, 

based storage scheme. In particular, the Figure shows a FIG. 1-18 of Appendix B is a diagram of the construction 

quadtree mesh with three levels of tree nodes and 10 leaf of an ADT for a data set with 5 points in a two dimensional 
nodes for the cells in the computational domain. 5 hvperspace. In particular, the Figure shows an ADT built for 

FIG. 1-9 of Appendix B are respective diagrams of data with two search keys. Each point in the data set is linked 

tree-traversal paths for locating the north face-neighbors of to the tree node which is used to store it. According to a 

particular cells. This provides information about a cell's well-known algorithm, a list of objects falling within a 

logical relationship to other cells in the domain. In specified target range is returned. According to one embodi- 

p articular, the FIG. 3 shows a cell-centered flow solver on a to ment of the present invention, the ADT is used for finding 
2-D mesh. Cell-centered storage of state data implies that intersection candidates between finite sized objects, 

each cell will exchange flux with its nearest face neighbors. Accordingly, a Cartesian bounding box of a finite sized 

The tree structure is relatively compact, and provides vir- object in d dimensions is defined as the smallest Cartesian 

tually all mesh connectivity information, while permitting region which can contain the object, 

direct computation of geometric properties. Geometry and 15 FIG. 1-19 of Appendix B are diagrams of Cartesian 
cell size need not be stored explicitly. Adaptation is rela- bounding boxes according to a number of predetermined 

lively easy to implement, since the structure is extensible embodiments for use according to the present invention. In 

ihiough division of leal nodes. Accordingly, it is possible to particular, the Figure show's examples for planar and 3-D 

produce a diagram of isobars and mesh cuts on a business jet objects. 

configuration computed with an octree based approach 20 FIG. 1-20 ol Appendix B is a diagram of inside, outside, 
accoiding to the prior art. In particular, it is possible to and intersected triangles in a predetermined domain cut by 

produce flow examples computed with octree and quadtree a boundary. In particular, the FIG. shows a sketch in which 

meshes. The mesh-cut through the aircraft symmetry' plane a set of triangles on an unstructured mesh are cut bv a closed 

provides an indication of mesh resolution. It is further boundary. It is desired to classify the triangles in the mesh 

possible to produce diagrams of density contours and 25 as inside or outside of the boundary. The^ intersected tri- 
adapted quadtree grids showing time histories of projectile angles have already been identified. A seed for a painting 

penetration, according to the prior art. algorithm is identified by determining the status of one of the 

FIG. 1-13 of Appendix B are diagrams of prior art data unclassified triangles. The seed triangle is classified by a full 

structures for polyhedral tessellations in two and three query evaluation and is painted according to the result, to 

dimensions. In particular, the Figure depicts common local 30 indicate its classification. Using this triangle as the starting 
data structures for describing mesh connectivity. One such point, the seed triangle paints each of its non-intersected and 

structure lists cell-to-vertex connectivity of an element, or unclassified neighbors with its own status. In turn, each of 

the cell-to-edge (i.e, cell-to-face) connections. The edges these triangles then recursively passes the result to its 

may he directed and may include links to cells on either side, non-intersected and unpainted neighbors, 
thus straddling mesh elements, and enabling communication 35 FIG. 2-2 of Appendix B is a diagram of surface triangu- 
therebetween. Further, winged edge structures are known lations for two intersecting teardrop shaped bodies, with the 

and can be used in connection with the present invention. labels indicating regions painted by a single seed. In 

Accordingly, it is possible to produce diagrams of adapted particular, the Figure shows two closed surface triangula- 

mesh, and computed isobars for in viscid flow over an tions which intersect tip-to-tail. 

ONERA M6 wing at particular predetermined conditions, 40 FIG. 2-1 of Appendix B is a diagram of component-based 
computed using an unstructured representation of the Car- Cartesian mesh generation, with and without removal of 
tesian mesh, and adapted Cartesian mesh and computed internal geometry. In particular, the Figure shows the pres- 

isobars for wing representations, as well as a diagram of an ence of internal geometry, causing intersection of a compo- 

adapted mesh and computed isobars for in viscid flow over a nent boundary by a body-cut Cartesian mesh cell. A cell may 

high wing transport configuration with an unstructured Car- 45 be completely encased within a component. Before being 
tesian mesh containing 2.9M cells with 10 adaptations. It is rejected from the mesh, the cell is tested for such contain- 

further possible to combine unstructured data storage with a ment. 

component-based surface modeling for complete aircraft FIG. 2-2 of Appendix B is a diagram of selected geomet- 
geometries. Surface triangulations are CAD generated and ric degeneracies in two and three dimensions. In particular, 

are permitted to overlap and intersect. It is possible to 50 the Figure shows degeneracies such as improper 
produce a diagram of a high wing transport configuration intersection, colinearity, and coplanarity, 

with a high-lift system deployed, and with the mesh con- FIG. 2-3 of Appendix B is a diagram of an intersecting 

taining 1 .65M cells at 10 levels of refinement and with the pair of generally positioned triangles in three dimensions. In 
mesh presented by cutting planes at 3 spanwise locations particular, the Figure shows a view of two intersecting 
and the cutting plane on the starboard wind being flooded by 55 triangles, with each intersecting triangle -triangle pair con- 
isobars of a discrete solution. tributing one segment to the final polyhedron which com- 

FIG. 1-17 of Appendix B is a diagram of the relationship prises the wetted surface of the configuration. With the task 

between partitions of a search space with partition directions of intersecting a particular triangle reduced to an intersection 

alternating over the keys of a search space, according to the test between that triangle and those on the list of candidates 

prior art. In particular, the Figure shows the partitioning of 60 provided by the ADT, the intersection problem is re-cast as 
a selected hyperspace into resulting sub-regions, to ensure a series of tri-tri intersection computations. The Figure 
that the children of each node contain the same amount of shows a view of two intersecting triangles as a model for 

objects. The tree accordingly remains balanced, even when discussion. Each intersecting tri-tri pair will contribute one 

the data is non-uniformly distributed in the search space. For segment to the final polyhedra that will comprise the wetted 

each node in the tree, the location of the partition associated 65 surface of the configuration. The assumption of data in 
with it is stored. The tree holds the maximum amount of data general (as opposed to arbitrary ) position implies that the 

at every level, and is therefore not as deep as an unbalanced intersection is always non-degenerate. Triangles may not 
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shore vertices, and edges of tri-tri pairs do not intersect 
exactly. Thus, all intersections will be proper. This restric- 
tion will be lifted in later sections with the introduction of an 
automatic tie -breaking algorithm. 

FIG. 2-4 of Appendix B is a diagram of the implemen- 
tation of a slope-pierce test on a selected edge against a 
triangle under test. In particular, the Figure checks for the 
intersection of a particular edge through a selected triangle. 
According to one approach, the intersection p of the segment 
with the triangle is evaluated for containment by verifying 
the same sign in each ot three cross-products. A conven- 
tional implementation of this approach has failures, as p falls 
near the edges or vertices of the triangle. 

FIG. 2-5 of Appendix B is a diagram of a Boolean test 
used in connection with the present invention to check 
whether a particular edge crosses a plane defined by a tested 
triangle, through computation of the signed volumes of two 
tetrahedra. In particular, the Figure shows that if a and b are 
established on opposite sides of the plane. It is thus deter- 
mined whether segment ab pierces within the boundary ol 
the triangle. This only occurs if the three tetrahedra formed 
by connecting the endpoints of ab with the endpoints of the 
three edges of the triangle 012 all have the same sign. 

FIG. 2-6 of Appendix B is a diagram of a Boolean test for 
the piercing of a particular line segment with the boundary 
of a triangle under test. In particular, the Figure illustrates a 
case in which all of the three volumes are positive. After 
determining the existence of all the segments which result 
from intersections between tri-tri pairs and connecting a 
linked list of all such segments to the triangles that intersect 
to produce them, all that remains is to actually compute the 
locations of the pierce points. This is accomplished by using 
a parametric representation of each intersected triangle and 
the edge which pierces it. 

FIG. 2-7 of Appendix B is a diagram of a parametric 
representation of a line defined by its endpoints and a plane 
defined by the vertices of a selected triangle in 3-D. Once a 
pierce point has been established to exist, a constructor is 
designed to generate the actual geometry of the pierce point. 

FIG. 2-8 of Appendix B is a diagram of respective 
configurations of intersecting cubes, and a wireframe show- 
ing intersection with segments linked to triangles. In 
particular, the Figure shows a configuration of two inter- 
secting cubes, as well as a wire-frame intersection of the 
cubes shown in highlighted form. 

FIG. 2-9 of Appendix B is a diagram of the decomposition 
of an intersected triangle using a constrained Delaunay 
triangulation algorithm, with the constraining segments 
shown as heavy solid lines. In particular, the Figure shows 
intersected triangles linked to 3 segments. The segments 
divide the intersected triangle into polygonal regions which 
are either completely inside or outside of the body. In order 
to remove the portions of these triangles which are inside, 
these polygonal regions are triangulated within each inter- 
sected triangle, permitting rejection of the triangles which 
lie completely inside the body. The final result of the 
intersection step is a list of segments linked to each inter- 
sected triangle. These segments divide the intersected tri- 
angles into polygonal regions which are either completely 
inside or outside of the body. In order to remove the portions 
of these triangles which are inside, we triangulate these 
polygonal regions within each intersected triangle and then 
reject the triangles which lie inside the body. The Figure 
shows a typical intersected triangle divided into two polygo- 
nal regions with the segments resulting from the intersection 
calculation highlighted. In the sketch, the two polygonal 
regions formed by the triangle's boundary and the segments 
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from the intersection have been decomposed into sets of 
triangles. Since the segments may cut the triangle arbitrarily, 
a pre-disposition for creating triangles with arbitrarily small 
angles exists. In an effort to keep the triangulations as well 
5 behaved as possible, we employ a two dimensional 
Delaunay algorithm within each original intersected tri- 
angle. Using the intersection segments as constraints, the 
algorithm runs within each intersected triangle, producing 
new- triangles which may be uniquely characterized as either 
inside or outside of the configuration. A variety of 
approaches to constructing the Delaunay triangulation ot a 
planar graph exist. However, since each triangulation to be 
constructed starts w ith the three vertices of the original 
intersected triangle, the incremental algorithm of Green and 
Sibson is appealing. Starting with the three vertices defining 
the original triangle, the pierce points are successively 
inserted into the evolving triangulation. After all the pierce 
points are inserted, the constraints are enforced using a 
simple recursive edge sw apping routine. 

FIG. 2-10 of Appendix B is a Voronoi diagram of a 
predetermined set with 35 sites. In particular, the Figure 
shows 35 sites having regions surrounding the sites identi- 
fying locations in a plane which are closer to that site than 
to any other site in a set. 

FIG. 2-11 of Appendix B is a diagram of a Voronoi 
triangulation of a set of 35 sites as shown in the immediately 
preceding Figure. Each triangulation to be constructed starts 
with the three vertices of the original intersected triangle. 

FIG. 2-12 of Appendix B is a diagram of a demonstration 
of a swap operation which exchanges the diagonal of the 
quadrilateral. In particular, the Figure shows identification 
of an edge to be swapped, and identification of the quadri- 
lateral formed by the two triangles that share the edge. The 
sw-ap operation reconfigures the quadrilateral by exchanging 
^ the diagonal. The Green and Sibson algorithm is extensively 
documented in the literature. Its central operation is the local 
application of an “incircle predicate" which determines 
which edges need swapping to make the triangulation 
Delaunay after each successive point insertion. This predi- 
4() cate examines the four points of a quadrilateral formed by 
two triangles which share a common edge. In the Figure, if 
the point d falls within the circle defined by (a,b,c) then the 
diagonal of the quad must be swapped (cb-»ad). Relating 
this discussion to the signed volume calculation starts by 
recognizing that if one projects the 2D coordinates (x,y ) of 
each point in the incircle test onto a unit paraboloid z=x 2 +y 2 
with the mapping, 

(k x . kyb-tik* k y k x 2 +k y 2 ) 

50 with k e {a,b,c,d} 

then the four points of the quadrilateral form a tetrahedron 
in 3D, and the incircle predicate may be viewed precisely as 
an evaluation of the volume of this tetrahedron. If V(T„\ 
b\c\d’):>0 then point d lies within the circle defined by 
55 (a,b,c) and the edge cb must be swapped to da for the 
triangulation to be Delaunay. 

FIG. 2-13 of Appendix B is a diagram of a Green and 
Sibson Delaunay triangulation algorithm. The algorithm 
begins with the intersected triangle which is to be retrian- 
60 gulated. A next site p is inserted and a triangle abc is located 
which contains the site. The new site is connected with the 
vertices of the triangle which contains it, by adding three 
edges pa, pb, and pc. These edges are incident upon p, and 
the original edges ac, ba, and cb are suspect. FIG. 2-14 of 
65 Appendix B is a diagram of incircle testing for a selected 
point z for containment within the circumcircle of p,q,r with 
z being contained and the diagonal of the quadrilateral pqzr 
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being swapped. In particular, the Figure shows the applica- segment-triangle intersection algorithm. The ADT returns 
tion of an incircle predicate to the suspect edges. If the test the list of intersection candidates while the logical tests 

point falls within the circle, the diagonal of the quadrilateral check for intersections. Counting the number of such inter- 

is swapped. Upon swapping, new suspect edges are identi- sections determines a triangle's status as inside or outside, 

fied on the boundary of the quadrilateral which are tested. 5 To avoid casting as many rays as there are triangles, a 
thus propagating the incircle test through the mesh. For “painting algorithm'' allows each tested triangle to pass on 

enforcement of constraining segments through edge deletion its status to the three triangles which share its edges. The 

and retriangulation, a particular constraining segment is in algorithm then recurses upon the neighboring triangles until 

the set ot edges ol the triangle. It it is not, then the edges it a constrained edge is encountered at which time it stops. In 

crosses are recursively swapped until the edge becomes an in this way the entire configuration is “painted" using very few 
edge in the triangulation. The edges for a segment that ray casts. The recursive algorithm is implemented using a 

intersects are removed and another is inserted. The vacant stack to avoid the associated overhead. Accordingly, it* is 

polygons on eithei side ol the segment are then retriangu- possible to produce a diagram of ray casting combined with 

lated. After constraint enforcement, each intersected triangle mesh painting with three internal regions shown identified 

is decomposed into small triangles which lie entirely inside is and painting with one ray per region, e.g. ray casting and 
oi outside of the wetted surface of the triangulation. painting with 3 teardrop bodies, Component triangulation 

Accordingly, it is possible to produce a diagram of decom- can be shown in wireframe, and three internal regions can be 

position of intersected triangles into triangles which lie painted. Each internal region can be seeded with a single rav 

completely inside or outside of the indicated configuration. cast, and filled in with the painting algorithm. 

FIG. 2-17 of Appendix B is a diagram of re triangulation 20 FIG. 2-24 of Appendix B is a diagram of unshaded 
within a large fuselage triangle pierced by a wing leading triangles using a 3x3 form of a simplex determinant without 

edge component with significantly higher resolution, the 52 incurring round-off error due to the initial subtraction of 

indicated segments describing the intersection of the leading coordinate data. In particular, the Figure shows an example 

edge constraining the triangulation. In particular, the Figure of surface triangulation on a teardrop shaped body. Triangles 

shows an example of re triangulation applied within a single 25 for which floating-point subtraction of their vertices is & not 
large fuselage triangle that has been intersected by a wing exact are shaded. To ensure accurate and robust computation 

leading edge. Since the wing leading edge has much higher according to the present invention, an orientation test is 

resolution requirements, its triangulation is substantially computed using data in floating point form. Further, the 

more refined. In all, 52 segments from the wing leading edge maximum round-off error bound is computed for the floating 

constrain the triangulation. Accordingly, it is possible to 30 point evaluation. Next, if the absolute magnitude of the 
produce a diagram of constrained Delaunay retriangulation signed volume from the orientation test is less than the error 

of intersected triangles for four intersecting teardrop bodies, bound, the orientation test is recomputed using exact, adap- 

e.g. a 3-D example of four intersecting teardrop bodies. The five precision floating point arithmetic, 

constrained Delaunay triangulations are shown with heavy FIG. 2-25 of Appendix B is a diagram of an expression 
lines - 35 tree for a 2-D orientation predicate. In particular, the Figure 

FIG. 2-18 of Appendix B is a diagram of point-in-polygon shows a list of approximate and true values in an expression 

testing using the (left) winding number and the (right) tree for evaluation of a determinant used to derive an error 

ray-casting approaches. In particular, FIG. 31 shows the bound for a 2-D form of the orientation test. Each true term 

testing of point q for containment within polygon P. A is expressed as an approximate term and an associated error 

winding number is computed by completely traversing the 40 term. Accordingly, the correctness of a sign determination is 
closed boundary P from the position of an observer at q, guaranteed mathematically. 

keeping a running total of signed angles between successive FIG. 2-26 of Appendix B is a diagram of an expression 
polygonal faces. Ray casting is simpler computationally. tree for three dimensional orientation, for computation of the 

Computation of a winding number involves floating point signed volume of a tetrahedron, Tabcd. In particular, the 

computation of small angles, which is prone to round off 45 Figure illustrates derivation of an error bound for this 
errors. The intersection and constrained triangulation rou- predicate. If the magnitude of the signed volume is less than 

tines of the two preceding sections have resulted in a set of a determined maximum value, then the predicate is deter- 

tri angles which may now be uniquely classified as either mined using exact arithmetic. The signed volume computa- 
intemal or exposed on the wetted surface of the configura- tion for arbitrarily positioned geometry can return a result 
tion. The only step left is then to delete those triangles which 50 which is positive (+1 ), negative (— 1 ) or zero (0), where +1 
are internal to the configuration. This is a specific applica- are non-degenerate cases and zero represents some geomet- 
tion of the classic “A point in polyhedron” problem for ric degeneracy. Implementation of this predicate, however, 
which we adopt a ray-casting approach. This method fits can be tricky, since it requires that we distinguish round-off 
particularly well within the framework of proximity testing error from an exact zero. Such considerations usually lead 
and primitive geometric computations. Simply stated, we 55 practitioners to implement the predicate with exact (arbitrary 
determine if a point p^p 0 ,p,,p 2 ) is within the ^component precision) arithmetic or with strictly integer math, 
of a configuration Q, if a ray, r, cast from p intersects the Unfortunately, while much hardware development has gone 
closed boundary C1 1 an odd number of times. The preceding into rapid floating point computation, few hardware archi- 
two sections demonstrated that both the intersection and tectures are optimized for either the arbitrary precision or 
triangulation algorithms could be based upon Boolean 60 integer math alternatives. In an effort to perform as much of 
operations checking the sign of the 3x3 determinant, and the the computation as possible on the floating-point hardware, 
same is true for the ray casting step. Assuming that r is cast we first compute in floating point, and then make an 
along a coordinate axis <+x for example) it may be truncated a-posteriori estimate of the maximum possible value of the 
just outside the +x face of the bounding box for the entire round-off error, e REmav . If this error is larger than the 
configuration T This ray may then be represented by a 65 computed signed volume, then the case is considered inde- 
line segment from the test point (p,Pi,p 2 ) to (ld^l v +e, p t ,p 2 ) terminate and we invoke the adaptive precision exact arith- 
and the problem reduces to a proximity query and the metic procedure developed by Shewchuk. If the case turns 
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out to be identically zero, we then resolve the degeneracy 
with a general tie-breaking algorithm based on a virtual 
perturbation approach. An error bound, €. K/ f , utv s is computed 
for floating point computation of the 3x3 determinant. The 
derivation accounts not only for the error in computing the 
determinant, but also for the error associated with floating 
point computation of the bound itself. This bound may be 
expressed as: 

?€+56€ >© < a ( $a w 0o. r i 

with: 
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consider the following two-dimensional version of the sim- 
plex determinant. 
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If the points a,b,c are assumed to be indexed with i=0, 1 . 2 
h) respectively, then taking 5=2 produces a perturbation matrix 
with: 
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where the circle ( ) overstrike on and x indicates that 
the operations are evaluated using floating point operations 
on IEEE 754 compliant hardware. € is precisely e=2 ; ' where 
p is the number of bits of the signifieand used by the 20 
machine, p may be evaluated by determining the largest 
exponent for which 1 .002"' =1 .0 when the sum and equality 
are tested w ith floating point. On most 32-bit platforms with 
exact rounding p=53 tor double precision and p=24 for 
single. In practice, only a very small fraction of the deter- 25 
minant evaluations fail. 

FIG. 2-27 of Appendix B is a diagram of real and 
perturbed configurations for degeneracy breaking, with per- 
turbations applied only to the vertical x 1 coordinate, for 
clarity. In particular, FIG. 35 shows an example of a per- 30 
turbation applied to a ray-casting determination in 2-D. A 
ray cast from point 5 intersected improperly wdth segments 
36, 47, 78, and 12, and is colinear with 64 and 81. By 
applying a perturbation according to the present invention, 
the indicated degeneracies are resolved as shown. 35 
Accordingly, point 6 is perceived as being below the line as 
point 6’, resulting in a tie-breaking effect. With degenerate 
geometry identified by the exact arithmetic routines, we can 
now remove the restriction imposed by the initial assump- 
tion that all input geometric data lie in general position. The 40 
richness of possible degeneracies in three dimensions cannot 
be overstated, and without some systematic method of 
identifying and coping with them, handling of special cases 
can permeate, or even dominate the design of a geometric 
algorithm. Simulation of Simplicity (SoS) is one of a cat- 45 
egory of general approaches to degenerate geometry known 
generically as “virtual perturbation algorithms”. The basic 
premise is to assume that all input data undergoes a unique, 
ordered perturbation such that all ties are broken (i.e. data in 
special position is perturbed into general position). When a 
tie is encountered, we rely on the perturbations to break the 
tie. Since the perturbations are both unique and constant, any 
tie in the input geometry will always be resolved in a 
topologically consistent manner. Since the perturbations are 
virtual, no given geometric data is ever altered. The pertur- 
bation e(ij) at any point is a function of the point's index, 

i €{0, 1 N— 1} and the coordinate direction,] e{l 

d}. According to the present invention, a perturbation of the 
following form is used: 

o < i < N - l 

£(/. j) = e where [ <: j <d 
8>d 

This choice indicates that the perturbation applied to i, is 
greater than that on k, iff (i<k) or (i=k) 6 (j>l ). To illustrate. 


Taking the determinant of the perturbed matrix M () =M+0 
yields: 

dct\M j = dvi\M | + V/v, + cn) + 

E ] : (/>i - < | ) + f- 'iiln ~ Oi) + 

tc ' : { 1 ) + | + ( 1 ) + 

1 } + higher order terms. 

Since the data, a,b,c span a finite region in 2-space, intu- 
itively one can always choose a perturbation small enough 
such that increasing powers of e always lead to terms with 
decreasing magnitude. If det[M] ever evaluates to an exact 
zero, the sign of the determinant will be determined by the 
sign of the next significant coefficient in the e expansion. If 
the next term also yields an exact zero, we continue check- 
ing the signs of the coefficients until a non-zero term 
appears. 

FIG. 2-28 of Appendix B is a diagram of two improperly 
intersecting right parallelepipeds with degeneracies resolved 
using virtual perturbations and exact arithmetic, respectively 
showing the components before intersection showing 
degeneracy, and the result after degeneracy and tie-breaking. 
In particular, the Figure shows large and small cubes abut- 
ting exactly. By tie breaking with virtual perturbations in 
accordance with the present invention, as well as computing 
the intersection, retriangulating, and extracting the wetted 
surface, the degeneracy is resolved. The virtual perturbation 
scheme according to the present invention, resolved not only 
the coplanar degeneracy, but all improper edge-edge 
intersections, resulting in a proper overlap of the two poly- 
hedra. The exact arithmetic and tie breaking routines accord- 
ing to the present invention provide a robust intersection 
algorithm including intersection, triangulation, and 
raycasting, even in the case of degeneracies. 

FIG. 3-3 of Appendix B is a diagram of a list of triangles 
associated with children of a cut-cell, which is obtained 
according to the present invention, using ADT, or by exhaus- 
tively searching over the parent cell s triangle list. In 
particular, FIG. 37 shows the passing of a parent cell’s 
triangle list to its children. The triangles describe the wetted 
surface of a selected configuration. When a cell is 
subdivided, a child cell inherits the triangle list of its parent. 
As the mesh subdivision continues, the triangle lists con- 
nected to a surface intersecting (i.e., cut) Cartesian cell will 
get shorter by approximately a factor of 4 with each suc- 
cessive subdivision. According to the present invention, it 
becomes successively computationally advantageous to con- 
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duct searches over a cell's triangle list for intersection rather 
than performing ADT. 

FIG. 3-4.</ and FIG. 3-4./? of Appendix B are a diagram of 
a measurement of the maximum angular variation within a 
cut-cell i. In particular, they show the determination of an 5 
internal angular variation value, based upon a surface nor- 
mal vector, within each cut cell. These additionally provide 
a diagram of a measurement of the angular variation 
between adjacent cut-cells. In particular, they show deter- 
mination of a cell-to-cell angular average variation value, to 
The respective internal and cell-to-cell angular variation 
values are used according to the present invention to deter- 
mine whether to tag cut cells tor geometric refinement. 
When the angles are identically zero, all the cut cells will be 
tagged for refinement. When they are 180 degrees, only is 
those cells at sharp cusps will be tagged. 

FIG. 3-5 of Appendix B is a diagram of the effect of 
varying angle threshold on number of cells produced after 
refinement of three different configurations. In particular, the 
Figure shows the sensitivity of the refinement to variation of 2(> 
parameters for angles from zero to 180 degrees, for three 
configurations. All surface intersecting Cartesian cells in the 
domain are initially automatically refined a specified number 
of times (R„„„) r By default this level is set to be 4 divisions 
less than the maximum allowable number of divisions 25 
<R„ im ) y in each direction. When a cut cell is tagged for 
division, the refinement is propagated several (usually 3-5) 
layers into the mesh using a “buffering" algorithm which 
operates by sweeps over the faces of the cells. Further 
refinement is based upon a curvature detection strategy. This 30 
is a two-pass strategy which first detects angular variation of 
the surface normal, n, within each cut cell and then examines 
the average surface normal behavior between two adjacent 
cut cells. Taking k as a running index to sweep over the set 
of triangles, T,, V ; is the j"' component of the vector 35 
subtraction between the maximum and minimum normal 
vectors in each Cartesian direction. 

\y=ma x(n L )\fik 

The min(-) and max(-) are performed over all elements of 40 
T,-. The angular variation within cell i is then simply the 
direction cosines of V : 

maxi*; ( ) - minK . ) 

cos 0, ) = 45 

r 

Similarly, (<j>y) r 5 measures the J rh component of the angular 
variation between any two adjacent cut cells r and s. With n,. 
denoting the unweighted unit normal vector within any cut 50 
cell i, the components of are: 


If 0j or in any cell exceeds an angle threshold (usually set 
to 25°) the offending cell is tagged for subdivision in 
direction j. 

FIG. 3-7 of Appendix B is a diagram of a Cartesian mesh 60 
with predetermined total divisions in each direction dis- 
cretizing using the region from xO to x l . In particular, the 
Figure illustrates the integer coordinate numbering scheme 
according to the present invention, in three dimensions. FIG. 

40 shows a model Cartesian mesh covering the region [x () , 65 
x, ]. Every cell in such a mesh can be uniquely located by the 
integer coordinates (i^ij.U) which correspond to the vertex 
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closest to x 0 . If we allocate m bits of memory to store each 
integer y, the upper bound on the permissible total number 
of vertices in each coordinate direction is 2"'. On a mesh 
with N ; prescribed nodes, performing R ; cell refinements in 
each direction will produce a mesh with a maximum integer 
coordinate of 2* / (N / — 1 H-l which must be resolved in m bits. 

2 /? '(A r i -l>+l £2'” 

Thus, the maximum number of cell subdivisions that may be 
performed in each direction is 

k )r4 log 2 ( 2'"- J >— log >( A r — I )J 

where the floor indicates rounding down to the next lower 
integer. Substituting back gives the total number of vertices 
which we can address in each coordinate direction. 

Af-2* i "‘ J '<A'-lRl 

Thus, M ; can never exceed 2"'. 

For any ceil i. is the integer position vector ( V {) ,V 0 ,V () ) 
If we also know the number of times that i has been divided 
in each direction, R ; ., we can express its other 7 vertices 
directly. 


F, = F 0 + < 0, 0.2*” 

*•'1 



F : = F 0 + ( 0, 2*" w *i 

-* i. 

0) 


F, = V„ + < o, ! ' 

•/?, 

Z Rnw 


V A = Fo + (2*' w 'o-*o. 

0, 

0) 


V s = F f) + 

0, 



F h = F 0 + {2 Rmax "~ R ", 


. 0) 


\/, - "> Rmti i j -R j 2-^2 ) 

Since the powers of two in this expression are simply a left 
shift of the bitwise representation of the integer subtraction 
Rmax^-R,, vertices V,-V 7 can be computed from V 0 and R y 
at very low cost. In addition, the total number of refinements 
in each direction will be a (relatively) small integer, thus its 
possible to pack all three components of R into a single 
32-bit word. 

FIG. 3-9 of Appendix B is a diagram of a general interior 
cell and first and second models for surface cut-cells. In 
particular, the Figure shows that at wall boundaries, Carte- 
sian cells are cut arbitrarily by the body geometry. The 
volume of each cut-cell which is inside the flow is computed 
by subtracting out the volume of the cell which protrudes 
into the wall. 

FIG. 3-10 of Appendix B is a diagram of a cut-cell in the 
abstract. In particular, the Figure shows the cut-cell linked to 
a set of four triangles which comprise the small swatch of 
the configuration's surface triangulation intersected by the 
cell. 

FIG. 3-11 of Appendix B is a diagram of an outcode and 
a facecode setup for a coordinate aligned region in two 
dimensions. In particular, FIG. 43 shows a two dimensional 
Cartesian cell covering a selected region. Particular points 
p-v are assumed vertices of the cells candidate triangle list. 
Each vertex is assigned an outcode associated with its 
location with respect to the cell. A central algorithm of any 
Cartesian mesh generation strategy involves testing the 
surface for intersection with the Cartesian cells. While the 
general edge-triangle intersection algorithm would provide 
one method of testing for such intersections, a more attrac- 
tive alternative comes from the literature on computer 
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graphics. This algorithm is highly specialized for use with 
coordinate aligned regions, and while it could be extended 
to non-Cartesian cells, or even other cell types, its speed and 
simplicity would be compromised. Since rapid cut-cell 
intersection is an important part of Cartesian mesh 5 
generation, we present a tew central operations of this 
algorithm in detail. More particularly, FIG. 43 shows a_two 
dimensional Cartesian cell c which covers the region [c,d]. 

The points (p, q v) are assumed to be vertices of c's 

candidate triangle list T . Each vertex is assigned an “out- |() 
code" associated with its location with respect to cell c. This 
boolean code has 2 bits for each coordinate direction. Since 
the region is coordinate aligned, a single inequality must be 
evaluated to set each bit in the outcode of the vertices. Using 
the operators & and | to denote bitwise applications of 
the □ and" and “or" boolean primatives. candidate edges 
(like rs) can be trivially rejected if 
outcode, & outcode A ^0 

Similarly, since outcode, |outcode v =0. the segment must be 
completely contained. If all the edges of a triangle, like tuv, , 0 
cannot be trivially rejected, then there is a possibility that it 
intersects the 0000 region. Such a polygon can be tested 
against the face-planes of the region by constructing a 
logical bounding box (using a bitwise “or") and testing 
against each facecode of the region. In FIG. 45, testing 25 
facecode, & (outcode, [outcode,, |outcode v ) 

Vje(0, 1, 2 2d- 1 ) 

results in a non-zero only for the 0100 face. When an 
intersection point, such as p’ or t' is computed, it can be 
classified and tested for containment on the boundary of [c, 30 
d] by examination of its outcode. However, since these 
points lie degenerately on the 01 XX boundary, the contents 
of this bit may not be trustworthy. For this reason, we mask 
out the questionable bit before examining the contents of 
these outcodes. Applying “not" in a bitwise manner yields: 35 
(outcode^, & rnfacecode t ) )= 0 , while 
(outcode,, & r^facecode,))^ 0 , 
which indicates that t’ is on the face, while p’ is not. 

FIG. 3-12 of Appendix B is a diagram of a divide-and- 
conquer strategy according to the prior art of Sutherland- 40 
Hodgman, with respect to polygon clipping, as a series of 
steps in which a polygon is clipped against a succession of 
infinite edges. In particular, the Figure shows a process for 
an arbitrary polygon clipped against a rectangular window 7 . 

FIG. 3-13 of Appendix B is a diagram of a setup for 45 
clipping a candidate triangle against a coordinate aligned 
region, including extraction of the clipped triangle. In par- 
ticular the Figure illustrates clipping for generating triangle- 
polygons. Accordingly, it is possible to produce a diagram of 
triangle-polygons on the surface of a high wing transport 50 
configuration resulting from the intersection of body-cut 
Cartesian cells with surface triangulation, with the triangle- 
polygons being triangulated and showing approximately 
500000 body-cut Cartesian cells. Accordingly, it is possible 
to produce a diagram of respective representations including 55 
cutting planes through a 4.72M cell Cartesian mesh, and a 
close-up of the mesh near the outboard nacelle, for a 
proposed supersonic transport design. Accordingly, it is 
possible to produce a diagram of the cutting planes through 
a mesh of multiple aircraft configurations with 5.6 1M cells 60 
and 683000 triangles in the triangulation of the wetted 
surface, showing portions of multiple cutting planes through 
a selected geometry, including for example a cutting plane 
at the tail of the rear two aircraft, or just under the helicopter 
geometry. 6? 

The detailed description of the preferred embodiments 
includes appendices A and B attached hereto and incorpo- 
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rated herein by reference in their entirety. More specifically, 
appendices A and B are as follows: 

Appendix A is M. J. Aftosmis, M. J. Berger, J. E. Melton, 
Robust ami Efficient Cartesian Mesh Generation for 
Component-Based Geometry, 35" J A1AA Aerospace 
Sciences Meeting and Exhibit. Jan. 6-9. 1997 (attached 
hereto and incorporated in the entirety. 

Appendix B is M. J. Aftosmis, Solution Adaptive Carte- 
sian Grid Methods for Aerodynamic Flows with Com- 
plex Geometries , Lecture notes for 28 r/l Computational 
Fluid Dynamics Lecture Series, von Karman Institute 
for Fluid Dynamics. Mar. 3-7, 1997 (attached hereto 
and incorporated herein in the entirety). 

What is claimed is: 

1. A computer implemented method for determining the 
intersection of selected component geometries where each 
component is described by a triangulated surface, eompiis- 
ing: 

locating triangle intersection candidates with a last spatial 
search: 

performing floating point filtering of geometric primitives 
by evaluating associated determinants in comparison 
with calculated error bounds, to determine whether an 
exact arithmetic determination of the determinants is 
necessary in view of the result of an evaluation being 
less than the error bound; 

resolving degeneracies by virtual perturbation; and 
identifying regions and minimizing ray casting by tra- 
versing connected regions. 

2. The method according to claim 1 wherein said alter- 
native digital tree data structure is stored in a selected 
computer memory. 

3. The method according to claim 1 wherein filtering 
includes making floating point error estimates to determine 
which floating point operations are to be computed exactly. 

4. The method according to claim 1 wherein adaptive 
precision software is used for unambiguous computation of 
geometric primitives. 

5. The method according to claim 1 wherein a degeneracy 
is resolved unambiguously. 

6 . The method according to claim 1 wherein a degeneracy 
is resolved using a simulation of simplicity technique 
according to one embodiment. 

7. The method according to claim 1 wherein regions are 
identified using painting. 

8 . The method according to claim 1 wherein ray casting 
is minimized using painting. 

9. A method of grid generation for Cartesian volume grids 
in three dimensions with embedded boundaries using a 
computer, comprising: 

performing a rapid triangle Cartesian cell intersection; 
employing an integer coordinate system in a multilevel 
Cartesian grid; 

using bitmasks for uniquely locating selected features of 
a cell; 

using outcodes for rapid intersection of segments of 
triangles and cells; 

matching split cells in adjacent polyhedra by testing 
adjacent Cartesian edge fragments for overlap; and 
identifying polyhedra in split cells. 

10. The method according to claim 9 wherein perfor- 
mance of a rapid triangle Cartesian cell intersection includes 
using an alternating digital tree data structure. 

11. The method according to claim 9 wherein said integer 
coordinate system is compressed. 
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12. The method according to claim 9 including constrain- 
ing selected features of a cell. 

13. The method according to claim 9 wherein said 
selected features are selected from a group of features 
including points, faces, edge, and corners. 

14. The method according to claim 9 wherein said polv- 
hedra are identified using painting. 

15. The method of determining intersection of compo- 
nents described by tri angulations in a selected configuration, 
comprising: 

performing a proximity search; 
conducting primitive geometric operations; 
performing adaptive precision exact arithmetic; and 
breaking geometric degeneracies. 

16. The method according to claim 15 wherein perfor- 
mance of a proximity search includes determining a list of 
triangles on components which may intersect with a triangle 
on a polyhedron under consideration. 

17. The method according to claim 16 wherein at least a 
single triangle from said list is compared for intersection 
with a triangle on a polyhedron under consideration. 

18. The method according to claim 17 including perform- 
ing signed volume tests on tetrahedrons constructed from 
said at least a single triangle and said triangle on a polyhe- 2 
dron. 

19. The method according to claim 15 including extract- 
ing the wetted surface of the intersected components. 

20. In both ray casting and Boolean tri angle -to- triangle 
intersection evaluations, the method comprising: 

performing floating point filtering; 
performing exact arithmetic; and 
performing virtual perturbation. 

21. The method according to claim 20 including perform- 
ing simulation of simplicity evaluation. 

22. The method according to claim 20 including evalu- 
ating the determinant of a matrix. 

23. A system for determining the intersection of selected 
component geometries where each component is described 
by a triangulated surface, comprising: 

means for locating triangle intersection candidates with 
fast spatial search; 

means for performing floating point filtering of geometric 
primitives by evaluating associated determinants in 
comparison with calculated error bounds, to determine 
whether an exact arithmetic determination of the deter- 
minants is necessary in view of the result of an evalu- 
ation being less than the error bound; 
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means for resolving degeneracies by virtual perturbation; 
and 

means for identifying regions and minimizing ray casting 
by traversing connected regions. 

' 24. The system according to claim 23 wherein said 

alternative digital tree data structure is stored in a selected 
computer memory. 

25. The system according to claim 23 wherein filtering 
includes making floating point error estimates to determine 

0 which floating point operations are to be computed exactly. 

26. The system according to claim 23 wherein adaptive 
precision software is used for unambiguous computation of 
geometric primitives. 

27. The system according to claim 23 wherein a degen- 

5 eracy is resolved unambiguously. 

28. The system according to claim 23 wherein a degen- 
eracy is resolved using a simulation of simplicity technique 
according to one embodiment. 

29. The system according to claim 23 wherein regions are 

* identified using painting. 

30. The system according to claim 23 wherein ray casting 
is minimized using painting. 

31. A system of generation for Cartesian volume grids in 
three dimensions with embedded boundaries, comprising: 

1 means for performing a rapid triangle Cartesian cell 

intersection; 

means for employing an integer coordinate system in a 
multilevel Cartesian grid; 

) means for using bitmasks for uniquely locating selected 
features of a cell; 

means for using outcodes for rapid intersection of seg- 
ments of triangles and cells; 

means for matching split cells in adjacent polyhedra by 
testing adjacent Cartesian edge fragments for overlap; 
and 

means for identifying polyhedra in split cells. 

32. The system according to claim 31 wherein perfor- 
mance of a rapid triangle Cartesian cell intersection includes 
using an alternating digital tree data structure. 

33. The system according to claim 31 wherein said integer 
coordinate system is compressed. 

34. The system according to claim 31 including constrain- 
ing selected features of a cell. 

35. The system according to claim 31 wherein said 
selected features are selected from a group of features 
including points, faces, edge, and comers. 

***** 


